% spnet.m: Spiking network with axonal conduction delays and STDP
% Created by Eugene M.Izhikevich.                February 3, 2004
M=100;                 % number of synapses per neuron
D=10;                  % maximal conduction delay 
% excitatory neurons   % inhibitory neurons      % total number 
Ne=800;                Ni=200;                   N=Ne+Ni;
a=[0.02*ones(Ne,1);    0.1*ones(Ni,1)];
d=[   8*ones(Ne,1);    2*ones(Ni,1)];
sm=5;                 % maximal synaptic strength
sm_1=1/sm;

load samples;

post=ceil([N*rand(Ne,M);Ne*rand(Ni,M)]); 
% for ii=1:400
% post(ii,201:400)=[402:2:800]';
% post(400+ii,201:400)=[2:2:400]';
% end;

dv=diff(samples(:,3));
dv(end+1)=dv(end);
dv(find(abs(dv)>0.2))=0;
 
s=[sm*rand(Ne,M);-sm*ones(Ni,M)];         % synaptic weights
%s=[6*ones(Ne,M);zeros(Ni,M)];
sd=zeros(N,M);                          % their derivatives
for i=1:N
  if i<=Ne
    for j=1:D
      delays{i,j}=M/D*(j-1)+(1:M/D);
    end;
  else
    delays{i,1}=1:M;
  end;
  pre{i}=find(post==i&s>0);             % pre excitatory neurons
  aux{i}=N*(D-1-ceil(ceil(pre{i}/N)/(M/D)))+1+mod(pre{i}-1,N);
end;
STDP = zeros(N,1001+D);
v = -65*ones(N,1);                      % initial values
u = 0.2.*v;                             % initial values
firings=[-D 0];                         % spike timings

for sec=1952:2200                      % simulation of 1 day
  for t=1:1000                          % simulation of 1 sec
    I=zeros(N,1);
    input=zeros(N,1);
    I(ceil (N*rand))=20;% random thalamic input 
    
    %if floor(t/10)==0
    if sec >100    
        % input torc
        %sm=10;
        %sm_1=1/sm;
        sind=rem((sec*100+floor(t/100)),10000);
        if sind==0
            sind=1;
        end;
        Iind=floor((samples(sind,5)/2/pi+1.1)*200);
        w=20;
        if Iind > w
            ind1=Iind-w;
        else
            ind1=Iind+1;
        end;

        for ii=ind1:Iind+w
            input(ii)=25*(1-sin((ii-Iind)/10/2/pi)^2);
        end;

        for ii=Iind+w:400
            input(ii)=1*rand(1);
        end;

        % output position
        %Iind=400+floor((samples(sind,3)/2/pi+1.1)*200);
        %dv(2:2000)=diff(samples(:,3));
        Iind=400+floor(dv(sind)*1000+200);
        w=20;

        if Iind -400 > w
            ind1=Iind-w;
        else
            ind1=Iind+1;
        end;

        for ii=ind1:Iind+w
            input(ii)=25*(1-sin((ii-Iind)/10/2/pi)^2);
        end;

        for ii=Iind+w:800
            input(ii)=1*rand(1);
        end;

        I(1:800)=input(1:800);

    end;

	if sec>2000
      I(401:800)=1*rand(400);
    end;
    
    fired = find(v>=30);                % indices of fired neurons
    v(fired)=-65;  
    u(fired)=u(fired)+d(fired);
    STDP(fired,t+D)=0.1;
    for k=1:length(fired)
      sd(pre{fired(k)})=sd(pre{fired(k)})+2.0*sm_1*STDP(N*t+aux{fired(k)}).*(sm-s(pre{fired(k)}));
    end;
    firings=[firings;t*ones(length(fired),1),fired];
    k=size(firings,1);
    while firings(k,1)>t-D
      del=delays{firings(k,2),t-firings(k,1)+1};
      ind = post(firings(k,2),del);
      I(ind)=I(ind)+s(firings(k,2), del)';
      sd(firings(k,2),del)=sd(firings(k,2),del)-1.5*sm_1*STDP(ind,t+D)'.*s(firings(k,2),del);
      k=k-1;
    end;
    v=v+0.5*((0.04*v+5).*v+140-u+I);    % for numerical 
    v=v+0.5*((0.04*v+5).*v+140-u+I);    % stability time 
    u=u+a.*(0.2*v-u);                   % step is 0.5 ms
    STDP(:,t+D+1)=0.95*STDP(:,t+D);     % tau = 20 ms
  end;
%figure(15);
%   plot(firings(:,1),firings(:,2),'.');
%   axis([0 1000 0 N]); drawnow;
%figure(16);
%hist(s(:),100);
  spikes=firings;
  disp([sec,length(firings(:,1)),min(min(s(1:800,:))),max(max(s(1:800,:)))]);

  if mod(sec,100)==0
   str=['sw' num2str(sec) '=s;'];
   eval(str);
   str=['save sw' num2str(sec) ' post sw' num2str(sec)]
   eval(str);
  end;

  if sec>1950 & sec<2200
str2=['sp' num2str(sec) '=spikes;'];
eval(str2);
str2=['save sp' num2str(sec) ' sp' num2str(sec)]
eval(str2);
end;

  STDP(:,1:D+1)=STDP(:,1001:1001+D);
  ind = find(firings(:,1) > 1001-D);
  firings=[-D 0;firings(ind,1)-1000,firings(ind,2)];
  s(1:Ne,:)=max(0,min(sm,0.01+s(1:Ne,:)+sd(1:Ne,:)));
  sd=0.9*sd;
%   cmat=zeros(1000,1000);
%   for ii=1:1000
%       cmat(ii,post(ii,:))=s(ii,:);
%   end;
%   figure
%   imagesc(cmat(1:800,:));
  
end;


